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Abstract. Fast Fourier transforms are used to develop algorithms for the fast generation 
of correlated Gaussian random fields on rectangular regions of R d . The complexities of the 
algorithms are derived, simulation results and error analysis are presented. 



1. Introduction 

In this note we present two algorithms for the fast simulation of Gaussian random fields 
(GRF) on a rectangular region of ~R d . 

Among the many possible applications we only want to mention here the simulation of 
pj^ numerical solutions of stochastic partial differential equations (SPDEs), and this was one 

^vq of the original motivations of the present note |10j . A typical example of an SPDE is the 

y—i stochastic heat equation 

d 

^ — u(t,x) = A x u(t,x) + rj(t,x), t > 0, x G D, 

in a domain D of M. d , driven by a centered Gaussian random field 7](t,x), which typically is 
"white", i.e., uncorrelated, in the time direction, and which has a covariance function C(x, y), 
x, y E D, in the spatial variable. Here A x stands for the Laplacian in the x variable. For 
accounts of the mathematical framework for SPDEs we refer the interested reader to, e.g., 
i— i El E2] and to the works cited there. 

In this paper we describe algorithms which generate approximations to centered Gaussian 
random fields on a rectangular grid with the help of the fast Fourier transform (FFT). This 
t**"*- entails that the resulting generated samples have periodic boundary conditions. If one uses 

instead discrete cosine and sine transforms (DCT and DST), similar results are obtained 
with Neumann and Dirichlet boundary conditions. A generalization to other domains than 
rectangles should be possible via NFFT (cf. (5j [Tlj Hj). An algorithm using FFT for fast 
Gaussian random field generation for the one-dimensional case can be found in [13]. Here 
t— I we generalize this algorithm to <i-dimensional rectangular regions. Furthermore we present a 

second algorithm that replaces one of the FFTs by a special data arrangement. In computer 
t> experiments this second algorithm is faster by a factor two. On an older computer the 

generation of one GRF of size 512 x 512 takes 0.25 seconds (250000 clocks). 

Spectral methods of different type can be found in [§] and references therein. Here we 
take the spectral representation of a given covariance, generate the Gaussian random field by 
multiplying the spectral density with the Fourier transformed white noise field and then apply 
the inverse Fourier transform. For an approach based on the circular embedding technique 
we refer the interested reader to [8] and the literature quoted there. 
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The paper is organized in the following way: Section [2] presents the properties of GRFs 
that will be used for the construction of the algorithms. These algorithms are presented in 
Section [3] and their complexity is analyzed. Finally Section [4] presents a class of covariance 
functions and shows simulation results and error estimates of an implementation in C++. 

2. Properties of continuous GRFs 

Assume that ip is a real-valued stationary Gaussian random field (GRF) on M. d . Thus 
for all x 6 W 1 , (p(x) defines a random variable, and the family (tp(x),x G M. d ) consists of 
identically distributed random variables on a probability space (f2, J-, P) that satisfy for any 
subset {xi, . . . , x n G and G K that the random variable Ylk=i ° k ^O^fc) * s normally 
distributed. Then its mean m = E(c^(x)) = E(y?(0)), x G and its covariance C(x) = 
K(ip(0)ip(x)) — m? = K(ip(y)ip(x + y)) — m 2 , x,y G M. d completely characterize the random 
field and are therefore its only relevant statistical parameters. Without loss of generality we 
will assume that m = 0. The goal of this section is to construct GRFs <f> for a given covariance 
C. We will present a heuristic approach here. A mathematically rigorous approach can be 
found in |10| . 

Let us start with the properties of the covariance function. We note that C is a positive 
semi-definite, symmetric function, and from now on we assume that it is continuous. Then 
Bochner's theorem [Tj states that C is the Fourier transform of a positive measure \xc on M. d , 
i.e. C can be written as 

C{x, y) = J e~ 2 ^-^ d^ c (p), y G R d , 

where (•,•) denotes the Euclidean inner product on M. d . We observe that C is actually a 
function of the difference x — y, and that it is an even function of x — y. For most practical 
purposes, there is no loss of generality, if we assume that nc has a Lebesgue density denoted 
by 7 which is even and positive. Then C can be written in the following way: 

C(x,y)= ( e- 27Ti( - p ' x - y) ~f(p) dp, x,y£R d . 
J«. d 

Let If be a Gaussian white noise random field on M. d also known as cylindrical Wiener 
process or Q— Wiener process with Q = 1, i.e. informally, W is a centered Gaussian family 
{W(x),x G M d } with covariance E(W (x)W (y)) = 5(x - y), x,y G R d . Set 

(1) <p(x) = (T- 1 j 1/2 TW)(x), x G R d , 

where T denotes the (f-dimensional Fourier transform and J 7-1 is its inverse. Then, since W 
is centered Gaussian, so is <p. The covariance of ip is given by 

E(<p(xMy)) 

= //// e ~ 2 ™ ((P '* )+( ^ )} ^) 1/2 ^^^^ dy'dpdq 

= I J " e ~ 2m{{P,X)Hq ' V) h(p) 1/2 j(q) 1/2 J e 2 ^ + ^'Ux'dpdq 

= j e- 27Ti ^ x - y h{p)dp 
= C(x,y). 
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Equation suggests our first algorithm to generate samples of 92 with given covariance C, 
and it is the first one presented in Section [3j 

Next we will work out a more efficient algorithm by replacing one Fourier transform with a 
faster operation. Our computer experiments with this algorithm took half the time of those 
done with the algorithm based on Equation ([I]). The goal is to construct TW directly, i.e. 
a complex GRF with the same distribution as J-W has to be generated. We consider the 
complex-valued random field J-W with 

TW{p) = I e 2ni(p ' x) W(x) dx 

for p G M. d . This is obviously centered Gaussian and the covariance is given by 

E{FW{p)FW{q)) = [[ e 27Ti « p ' x) ~ (q ' y)) E(W(x)W{y)) dxdy= [ e 2 ^~^ dx = S(p - q), 



and similarly 

E{FW{p)FW{q)) = 5(p + q). 
The following lemma gives a direct construction of J-W . 
Lemma 1. Let V be the complex-valued random field defined by 

YieV = -K + W and ImV = ir~W, 

where 

K + W{p) = I (W(p) + W{-p)) and TT'Wip) = \ (W{p) - W{-p)). 
Then V and TW have the same law. 

The lemma is proved by the following observations. As both random fields are centered 
Gaussian, one only has to show that they have the same covariance. First we calculate 

E{7r + W{p)7r-W(q)) 

= \ E((W(p) + W(-p))(W(q) - W{-q))) 

= \ {E{W{p)W{q)) + E(W{-p)W{q)) - E(W(p)W{-q)) - E(W{-p)W{-q))) 

= I (<KP " V) + S(p + q)- 5(p + q)- 5{p - q)) = 0. 

Therefore tt + W and tt~W are uncorrelated, and hence independent. Then this observation 
yields 

E(V(p)V(q)) = E(TT + W(p)7r + W(q)) - E(ir-W{p)Tr-W(q)) 

= IW(P ~q) + 25(p + q)- (25{p -q)- 26(p + q))) = 5{p + q) 

and similarly 



E(V(p)V(q))=5(p-q), 

and V = J-W in distribution as claimed in Lemma [TJ □ 
By Lemma [T] a random field tp given by 

(2) ip(x) := (J r-1 7 1 / 2 (7r + Ty + m~W)){x) 

is equal in distribution to ip, and in particular it has covariance C, too. So Equation ([2]) 
yields a second algorithm for the construction of random fields with covariance C which is 
also presented in Section [3j For the following we remark that V as defined above satisfies 
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V{—x) = V[x) for all x G Furthermore, by the previous calculations, Re V(x) and 

ImV(x) are independent centered Gaussian random variables. 

In the remainder of this section we shall consider the problem of discretizing the random 
fields above. Let A be a rectangular region in IR rf defined as the Cartesian product of d 
closed intervals Jj = [oi,6j], i = 1, d. For each i = 1, d, choose N G N points 
Uq < y\ < . . . < 2/7^-1 m \ a i-> bi) defining a partition of Jj. From now on fix N = (Aq, . . . , A^), 
and suppose that every N is even. We denote by K the set of all vectors K = (k\, . . . , kd) 
with h% = 0, . . . , Ni — 1. Define xk to be the point in A whose i— th coordinate (xx)i is equal 
to y^... Thus the xk, K G /C^, form a rectangular grid in A For K G /C^, define the cell A^- 
in j4 by 

d 

Ax = X [x K ,X K +eJ, 
i=l 

where &i = = 1, ...,d). Thus the set {A^, K G /C^} of cells forms a partition of 

A' = [ai,6i) x ... x [a d ,b d ). 

As above, we let W denote a white noise random field on IR d , which defines a discretized 
white noise random field W N by 

W N = (Wk , K G JC N ) 

with 

= lA^r 1 f W(x)dx. 
Ja k 

Hence the family {W^} is an independent family of centered Gaussian random variables, 
Wj[ having variance |A^| _1 . We may also view W N as a discrete white noise random field 
on A' by considering it as being constantly equal to on the cell Ak- 

Now we impose periodic boundary conditions at the boundary of A as will be automatically 
implemented by use of the FFT in our algorithms in the next section. That means, in the 
sequel we work on the torus defined by A. Therefore all calculations will be done modulo 
the side lengths of the rectangle A, or modulo the vector N, respectively. In particular, the 
following discrete analogues of Re V, Im V in Lemma [l] are well-defined for all K G K, N 

R N K 4« + w» K ), If = \ « - w? K ). 

Set V N = (Vg, K G JC N ) with 

t/JV _ r>N , ■ T N 

Observe that we have = — for all K G JC N , and therefore 
(3) V? K = Vg. 

In particular, we find that for all K G JC which are such that K = —K (mod N), Vg is 
real. Figure [l] shows an equidistant grid with = (4, 4) on the unit square of R 2 . The solid 
points are the grid points, and those which are blue have the property K = —K (mod (4, 4)). 
Consequently the random field is real at the blue grid points. The red points form an example 
of two points associated via K and —K with each other: K = (2,1) = —(2,3) (mod (4,4)). 
Thus the value of the random field at one of the red points is equal to the complex conjugate 
at the other. 
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Figure 1. A 4 x 4 grid on the unit square in 



For K £ JC N with K / — K, the independence of W£ and Wf! K entails 
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On the other hand for K £ JC N with K 



-K we get 
I A* 



E({R»f) - 1 ^ xr ' 



and — as already mentioned above — 1^ =0. 

It is straightforward to check that for every K S K N there is exactly one vector L £ K, 
so that L = -K (mod iV). Let £ 



N 



N 



denote any maximal subset of /C , so that the random 
variables Vff , X £ are independent. Thus we only have to generate the random variables 



and 1^- for K £ C N , and the values of the random field V at the other indices are 
determined from these by complex conjugation. In the appendix it is shown how one can 
construct one such set C N for an arbitrary dimension d £ N. Here we only remark that for 
d = 2, a possible choice for C N , iV = (N±,N2), is given by 

£ W,Ay ={k 1 = 0,..,,M 1 ,k 2 = 0,...,M 2 } 

a {fei = 1, . . . , Mi - 1, k 2 = M 2 + 1, . . . , 2M 2 - 1}, 

where we have set Mj = 2Vj/2, i = 1, 2. In Figure [2] this set is shown for an equipartition of 
the unit square with an 8 x 8 grid as those corresponding to the magenta points. The values 
of the random field at the green and black points is determined from those at the magenta 
points via Relation 

For simplicity we assume from now on equidistant partitions in each direction, and let 
| A^l denote the volume of the cells. Then the calculations above can be summarized by the 
following rules: 
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Figure 2. C N for an 8 x 8 grid given by the magenta points 



a) In the cell A^, K G /C, the value of the discrete random field to be simulated is given 
by 

(4) DFT- 1 7 (PK) 1/2 K + ^)|A|- 1 , 

where DFT denotes the discrete Fourier transform. 

b) the random variables R^, 1^, K G C N , are independent of each other; 

c) if K G C N is such that K = -K (mod N), then R% has the law Af(0, \ A N \- 1 ), 
T N - 0- 

d) if K G is such that K / —if (mod N), then i?^ and 1^ are distributed with the 
law Af(0, (2| A^l)- 1 ), and I? K = -1%; 

e) for K G IC, the points pk are chosen such that {px)i = i • {h — ai), for i = 1, 2. 

The corresponding algorithm implementing these rules in an efficient way can be found in the 
next section. 

We close this section with the following 

Remark 2. For higher dimensions d than d = 2, the implementation of C N becomes a rather 
tedious task, see the construction of C N in the appendix. Consider again the case d = 2 in 
figure^ If we add to C N the set of green points, i.e., some of the points at the boundary of the 
upper left subcube, then this set can be enumerated in a very simple way as k\ = 0, . . . , M±, 
&2 = 0, . . . , 2M2 — 1. After generation of the random variables, the values at the green points 
have to be overwritten according to Relations ([3]). Since the random variables which are 
superfluously generated are those at the boundary of the upper left cube, their number is 
negligible as compared to those in the interior of that cube. A similar consideration holds 
in the case of general d G N: If one generates additionally to the random variables indexed 
by C N those at the boundaries of the subcubes, one obtains a simple code at negligibly 
higher computing costs. After generation of all values, those which have been superfluously 
generated have to be overwritten with the correct values according to Equation ([3]). 
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3. Algorithms 

The following algorithms are based on the calculations of the previous section. They allow 
a fast generation of d-dimensional stationary Gaussian random fields with given covariance 
using the advantage of fast Fourier transforms instead of directly calculating computational 
expensive convolutions. The boundary conditions implemented here are periodic. For Neu- 
mann or Dirichlet boundary conditions, fast cosine and sine transformations have to be used. 

The first algorithm is an implementation of Equation Q. 

Algorithm 3. 

Remarks: 

(1) The functions FFT and FFT -1 include all necessary rescaling depending on the used 
FFT algorithm and the integers iVj. 

(2) A is a d-dimensional complex-valued array, B is real-valued. 

(3) x/ Cl ...fc d denotes the grid point corresponding to the integers (k±, . . . , kd)- The grid 
points are distributed equidistantly in each direction, i.e. the distance of two arbitrary 
neighbor grid points in direction e, is given by a constant Axj. 

(4) The points Pk v -k d i n the Fourier domain are given by (pki-k d )i = (^i ~~ Aj/2)//j for 
i = 1, . . . , d. 

Input: 

(1) d-dimensional rectangular region D, where 1%, . . . , Id is the length of the edges, 

(2) N\, . . . , Nd number of discretization points in each direction, all even, 

(3) 7 1 / 2 a symmetric, positive function on ~R d , 

(4) R() a function that generates independent AA(0, | A Ar | _1 )-distributed random numbers. 
Output: GRF B on D, where the covariance is given by the Fourier transform of 7. 

for ki = 0, . . . , iVj — 1, i = 1, . . . , d do 

B(ki, ...,kd)<— R(); 
end for 
A <r- FFT B; 

for ki = 0, . . . , iVj — 1, i = 1, . . . , d do 

A(ki, ...,k d )^ A{k x , ..-,k d )- 7(Pfci-fcJ 1/2 ■ PI" 1 ; 
end for 
B <- FFT" 1 ^; 

The second algorithm is based on Equation Q with Remark [2j For the direct construction 
of the Fourier transform of the white noise field, a random field with even real and odd 
imaginary part has to be constructed. At all grid points indexed by K G C N with K = —K 
mod(iV) the random field is real-valued as the random field is equal to its complex conjugate 
at these points. At the grid points where the random field is real-valued it is scaled with a 
factor one instead of 2" 1 / 2 . If we use a d-dimensional DFT on a grid with JVj discrete points 
in direction e«, the following algorithm generates a GRF with approximate covariance C. In 
order to keep the error small, one has to take care that the correlation length defined by the 
covariance is small in comparison to the size of the rectangular region. 

Algorithm 4. 

Remarks: 

(1) All calculations have to be done modulo Ni in the i-th direction. 
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(2) The function FFT includes all necessary rescaling depending on the used FFT 
algorithm and the integers TVj. 

(3) A is a d-dimensional complex-valued array, B is real-valued. 

(4) Xk v -k d denotes the grid point corresponding to the integers (fei, . . . , k d ). The grid 
points are distributed equidistantly in each direction, i.e. the distance of two arbitrary 
neighboring grid points in direction is given by a constant Axi. 

(5) The points Pki-k d hi the Fourier domain are given by (pki-k d )i = fti — Ni/2)/li for 
i = 1, . . . , d. 

Input: 

(1) d-dimensional rectangular region D with l±, . . . , l d length of the edges, 

(2) N± , . . . , N d number of discretization points in each direction, all even, 

(3) 7 1 / 2 a symmetric, positive function on R d , 

(4) RQ a function that generates independent J\f(0, (A^j "^-distributed random numbers. 
Output: GRF B on D, where the covariance is given by the Fourier transform of 7. 

for ki = 0, . . . , Nt - 1, i = 1, . . . , d - 1, k d = 0, . . . , N d /2 do 
if h G {0, Ni/2}, for all i = 1, . . . , d then 

ReA(k u ...,k d )^RQ- 7 (p fcl ... fc J 1/2 • ID]' 1 ; 

Im4(fci, . . . , k d ) <- 0; 
else 

ReA(k ± , ...,k d )^ 2-V 2 RQ ■ j{p kl -k d ) 1/2 ■ l^" 1 ; 

ImA(fei, ...,k d )<- 2- 1 / 2 R() ■ 7(p kl -k d ) 1/2 ■ PI" 1 ; 

ReA(N 1 -k 1 ,...,N d -k d )^ReA(k 1 ,...,k d ); 

ImA(JVi -k!,...,N d -k d )<- -lmA(k u k d ); 
end if 
end for 
B «- FFT^ 1 ^; 

Next, we give an algorithm that generates just the necessary amount of random numbers. 
For that algorithm, we observe that the two hyperplanes induced by ki = 0, . . . , JVj — 1 for 
i = 1, . . . , d — 1 and k d = or k d = N d /2 have the same structure as the corresponding torus 
in d — 1 dimensions. So we use a recursive implementation that solves the problem of setting 
the values of the random field at the grid points by setting all values except those in the two 
hyperplanes, and then restarts the algorithm for these hyperplanes in d — 1 dimensions. 

Algorithm 5. 

Remarks: 

(1) All calculations have to be done modulo Ni in the i-th direction. 

(2) The function FFT -1 includes all necessary rescaling depending on the used FFT 
algorithm and the integers 7Vj. 

(3) A is a d-dimensional complex-valued array, B is real-valued. 

(4) Xk l ...k d denotes the grid point corresponding to the integers (k±, . . . , k d ). The grid 
points are distributed equidistantly in each direction, i.e. the distance of two arbitrary 
neighbor grid points in direction ej is given by a constant Axj. 

(5) The points Pk\-k A h 1 the Fourier domain are given by (pfc 1 ...fc d )i = (ki — iVj/2)/Zj for 
i = 1, . . . , d. 

Input: 
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(1) (i-dimensional rectangular region D with li, . . . ,l d length of the edges, 

(2) N\, . . . , N d number of discretization points in each direction, all even, 

(3) 7 1 / 2 a symmetric, positive function on 

(4) R() a function that generates independent Af(0, | A Ar | _1 )-distributed random numbers. 

Output: GRF B on D, where the covariance is given by the Fourier transform of 7. 

function create -random _field(D , N±, . . . , N d ) 
A <— complex_array(-/Vi, . . . , N d ); 
set_array({N 1 ,...,N d }, { }, A); 
B <- FFT" 1 ^; 
end function 

function set-array ({Ni, . . . , Nj}, {kj + \, . . . , k d }, A) 
for ki = 0, . . . , N i} i = 1, . . . , j - 1, kj = 1, . . . , Nj/2 - 1 do 
ReAih, . . . , k d ) <r- 2- 1 /2 jR () . 7 ( Pfci ... fc Ji/2 . \ D \-l. 

ImA(%, ...,kd)4r- 2- 1 ' 2 R() ■ j( Pkl -k d ) 1/2 ■ \D\~ l ; 
Re^(iVi -k l ,...,N d -k d ) ^ ReA(h, ...,k d ); 
JmAlN 1 -ki,...,N d -k d ) -Im A(k x , . . . , k d ); 
if (\{N 1 ,...,N j }\ = l) then 

Re^(0, ...,kd)<- RQ ■ 7(P0fc 2 -fcJ 1/2 • l^l" 1 ; 
hnA(0,...,k d )<-0; 

Re A(N x /2, ...,kd)<-RQ- l(p Nl /2k 2 -k d ) 1/2 ■ PI -1 ; 
Im^iVxA...,^) ^0; 
else 

set-array{{Ni, iVj-i}, {0, fc^+i, . . . , A); 

seLarray({Ni, iVj-i}, {^/2, fc^+i, . . . , k d }, A); 
end if 
end for 
end function 

Finally we analyze the computational costs of the algorithms in this section. Let N = 
Ni + • • ■ + N d . Algorithm [3] sets the random array in O(N), each FFT is done in 0(iVlog N) 
and the costs for the multiplication with 7V2 is each for the real and the imaginary part O(N). 
So overall we have 0(iVlog N). The second algorithm [4] performs the generation of the com- 
plex random array in O(N) for real and imaginary part and then one FFT in O(NlogN). 
The additional generated random numbers cost 0(Ni + • • • + N d —i). Therefore the algorithm 
also has complexity 0(iVlog./V). Nevertheless, the constants are smaller. In simulations, 
Algorithm [3] took twice as long as Algorithm |4j Finally, Algorithm [5] constructs the com- 
plex random field in O(N) and performs the FFT in O^logA^), which overall leads to a 
complexity of 0(AnogAQ. Therefore all algorithms are in the same complexity class but 
in simulations they perform with different speeds due to different constants. So in terms of 
complexity, there is no gain if one omits one FFT, but in practice a factor two in speed could 
be of advantage. 

The algorithms have been implemented and tested for d = 1,2. The results are presented 
in the subsequent section. 
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(a) m = W,n = l,k = 1,1 = 2 



(b) m = 0.75, n = 1, fe = 1,1 = 2 





(c) m = 10, n = 3, fc = 1, i = 1 (d) m = 5, n = 1, A = 3, / = 1 

Figure 3. Simulation results with different covariance functions. 



4. Simulations and Statistical Tests 

This section shows the results of an implementation in C++ of the algorithms presented 
in Section [3] using FFTW [6j. The simulations were done on a rectangle in M? and on the 
interval [—it, it] in M. We first consider the two-dimensional implementation. Afterwards we 
do error analysis on the interval. It turned out that both algorithms lead to the same results 
but that Algorithm [4] is twice as fast as Algorithm [3} 

For the two-dimensional implementation, the covariance was chosen to be 



C(x,y) 



-ikUpp-v)^ dpj x,y £ 



with 



[m 



kl +(pf + P f) l r 



(5) 7 (p) 

where k, I, n G N, m G M>o, and p = (pi, P2) £ For k = I = n = 1 this covariance function 
is used in Euclidean quantum field theory, cf. [7]. It is shown in |10j that C has exponential 
decay for \x — y\ 3> 1. (This is well-known for k = I = n = 1, e.g. |7].) Statistical tests were 
made on the samples by choosing one grid point xq and calculating the statistical covariance 
of the corresponding random field with respect to the random field at all other grid points. 
Some results are shown in Figure [3| We remark that the graph in Figure 3(c) shows negative 
covariance values, however this is the correct behaviour and not an artifact of the simulation. 
A physical interpretation of that phenomenon would be some kind of countermove against the 
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displacement of a point like an elastic band. Figure |3(d)| shows a covariance function that is 
not symmetric under rotation. Functions like that might be interesting for applications with 
direction dependent correlations. 

The simulated random fields are visualized in Figure |4j [5j and [6| Two different possibilities 
how to visualize the data are presented in Figure |4| On the left, the random numbers are 
plotted as the graph of a function from R 2 to R. The figure on the right-hand side presents 
the random numbers in the form of colors where blue stands for small numbers and red 
for the large ones. Figure [5] uses this hot color map to compare correlated random fields 
that use the same white noise samples for W but different covariance functions. There are 
differences between all of the pictures, but it attracts attention that there are hardly any 
differences between the three pictures where 7 is of degree —4, while Picture |5(a) where 7 



is of degree —2 is completely different which is due to infinite variance of the underlying 
random field. This phenomenon can be observed for all tested other degrees as well. Figure [6] 
shows that using the same white noise and the same type of covariance function with different 
degrees of 7 leads to similar samples of the random fields but the higher the degree of the 
polynomial, the smoother the resulting noise as expected. Setting m to something larger 
than 1 will result in faster decreasing correlations which can be seen in Figure [3} 

We analyze the convergence of the statistical covariance on the interval [— ir, it] C R with 
theoretical covariance C defined by 

, , 32 • 10 6 1 

1\P) = 



vr (100 + p 2 ) 4 ' 

for pGK. The constant is used to scale the variance to 1. Then C can be computed explicitly, 
and it is given by 

C(x, y) = (200/3|x - y| 3 + 40|x - y\ 2 + 10|z - y\ + 1) exp(-10|x - y\), 

for x,y G R. For the error analysis we simulated 10 8 samples (rjj,j = 1,...,10 8 ) on the 
equidistant grids D n = {x? = i ■ vr/2 n - 1 ,i = -2 n -\ . . . , 2 n ~ 1 } for n = 2, . . . , 12. The 
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(a) 7 (P) = (1 + M 2 r 



(b) 7 (p) = (i + N 2 )- 





(c) 7 (p) = (1 + |p| 4 )- 



(d) 7(P) = (1+Pi+P2)" 1 



Figure 5. Images of size 512 x 512 with the same white noise. 






(a) 7 (p) = (1 + |p| 2 )- 



(b) 7 (p) = (1 + |p| 2 )- 



(c) 7 (p) = (1 + |p| 2 )- 






(d) 7 (p) = (1 + |p| 2 )- 



(e) 7 (P) = (1 + NT 



(f) 7 (p) = (1 + |p| 2 )- 



Figure 6. Images of size 512 x 512 with the same white noise. 
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-3 -2 -1 



space interval 

(a) theoretical covariance 




-3 -2 -1 



(b) 9 grid points 





(c) 33 grid points 



(d) 4097 grid points 



Figure 7. Theoretical and statistical covariance on different grids. 



covariance was then estimated by 



Cov(x r : 



10 s 

i=i 



for all i = — 2 ra_1 , . . . , 2 n_1 and n = 2, . . . , 12. Figure [7] shows the exact covariance C and the 
statistical results on a grid with 9, 33, and 4097 points. 

The error was computed by taking on each grid the maximum over all grid points of the 
difference of the theoretical and statistical covariance, i.e. the error e n on D n was computed 
by 



max 

i=-2"- 1 ,...,2 ri 



\Cov(x?)-C(x?)\. 



The results for n = 2, . . . ,6 are shown in Figure [8j The reference slope shows that the order 
of convergence is at least of order 0(iV -4 ) where N = \D n \. Finer grids are excluded from 
the plot because then the error is dominated by the error of the Monte Carlo sampling. 
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number of grid points 



Figure 8. Maximal difference of the statistical and theoretical covariance 
over all grid points and the reference slope of 0(iV~ 4 ). 



Appendix A. Construction of C n 

We suppose that d G N, and that TV = (N ± , . . . , N d ) with even N 1: N d G N. Put 
Mi = Ni/2, i = l, d. Recall that 

K N = {(k u ...,k d ), k i = 0,...,N i -l,i = l,...,d}, 

is subject to addition mod(iV). Set 

Co = {k G JC N , h = or Mi for alU = 1, . . . , d}. 

For n = 1, . . . , d construct as follows: Take all partitions (ji, . . . , j n ), (l n +i, . . . ,l d ) of 
(l,...,d), ordered in the natural way, and collect k G K N such that 

k n = l,...,M n -l, 

k n = 1, . . . , Mj. - 1, M j% + l,...,Nj.-l, i = l,...,n, 
k h =0,M h , i = n + l,...,d. 

Then set 

c N = y ^. 

n=0 

Observe that in the enumeration above, each /c G £ w appears precisely once. Moreover, the 
subset )Cq of K N so that —JCq = )Cq coincides with the subset Cq of C N . Set 



C N = C N \C$= |+J 



n=l 
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It is not hard to see that if k € C N , then —k £ K but — A; ^ jC^. Moreover, a simple 
counting argument shows that C N W (—C N ) has the same number of elements as /C^. Thus 
we find that K N = C$ l±J £ w 1+) (—£ N ), and hence C N has indeed all required properties. 
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